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Abstract 

We consider a PDE system comprising compressible hydrodynamics, flux-limited 
diffusion radiation transport and chemical ionization kinetics in a cosmologically- 
expanding universe. Under an operator-split framework, the cosmological hy- 
drodynamics equations are solved through the Piecewise Parabolic Method, as 
implemented in the Enzo community hydrodynamics code. The remainder of 
the model, including radiation transport, chemical ionization kinetics, and gas 
energy feedback, form a stiff coupled PDE system, which we solve using a fully- 
implicit inexact Newton approach, and which forms the crux of this paper. The 
inner linear Newton systems are solved using a Schur complement formulation, 
and employ a multigrid-preconditioned conjugate gradient solver for the inner 
Schur systems. We describe this approach and provide results on a suite of test 
problems, demonstrating its accuracy, robustness, and scalability to very large 
problems. 
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Methods, Implicit Methods 
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1. Introduction 

A fundamental physics component in cosmology and astrophysics applica- 
tions is the transport of ionizing radiation along with the interaction of that 
radiation with the hydrodynamic motion and ionization state of the surround- 
ing gas. One example currently receiving a great deal of attention is cosmic 
reionization our motivation in this work. Observations indicate that an 
early population of UV emitting galaxies photoionized the intergalactic hydro- 
gen and helium gas when the universe was about 800 million years old (rcdshift 
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~ 8). A computational challenge is to calculate this process self-consistently, 
coupling the radiative transfer of ionizing photons, the ionization kinetics and 
photo-heating of the gas, and the attendant hydrodynamic motions. This prob- 
lem is challenging because the physics is numerically stiff and cosmic reionization 
is intrinsically three-dimensional, involving the growth, percolation, and overlap 
of ionization zones around an irregular and evolving distribution of galaxies with 
time-dependent luminosities. In addition, the problem inherits the large range 
of spatial scales (~ 10 5 ) intrinsic in galaxy formation simulations, necessitat- 
ing the use of spatially adaptive mesh or particle-tree methods and large-scale 
parallel computing 0, Q . 

A variety of 3D radiative transfer methods are under development to tackle 
this problem 0]. These necessarily simplify the description of the radiation 
field in order to render the problem tractable. These methods include ray trac- 
ing using long and short characteristics [!, @, 0, 0] , Monte Carlo 0, [Tc| , and 



moment methods [11|, [12J, [13|. However, only some of these codes allow so- 



lution of the coupled problem on spatially-adaptive grids, and very few allow 
distributcd-mcmory parallelism. More importantly, all of these codes handle 
the interactions between hydrodynamic, radiative, and chemical processes in 
an explicit, operator-split fashion, which ignores stiff couplings that often arise 
between these components. When this happens, such codes must dramatically 
reduce allowable timesteps, or resort to subcycling, to ensure stability and ac- 
curacy of the coupled simulations. 

Radiation transport and chemical ionization time scales are much faster 
than typical hydrodynamic time scales. This is particularly evident in dense 
gas bound to galaxies where recombination and light crossing times are short 
compared to the age of the universe (Hubble time). Moreover, these processes 
are very tightly coupled since radiation induces ionization that in turn affects 
opacities. While time-explicit, adaptive subcycling schemes have been developed 
that are capable of returning accurate solutions in all regimes of interest Q , it is 
our view that for both computational efficiency and solution accuracy, tightly- 
coupled implicit methods require investigation. Here we present such a method. 

We solve ionizing radiation transport, chemical ionization kinetics, and gas 
photo-heating using a fully implicit inexact Newton method. This algorithm is 
coupled to a cosmological hydrodynamics code through an explicit, operator- 
split formalism. The inner linear Newton systems are solved using a Schur 
complement formulation, which neatly decouples the local microphysics from 
the transport calculation. Radiation transport is modeled in the flux-limited 
diffusion approximation for simplicity, although our approach can be easily gen- 
eralized to higher-order moment schemes. The use of Schur complements allows 
the application of optimal and scalable multigrid methods for the solution of 
the scalar radiation diffusion equation. We describe our algorithm in detail. 
We then illustrate our method's accuracy, robustness, and parallel scalability 
against a suite of verification tests of increasing size and complexity. In its 
current implementation, we are restricted to calculations on uniform Cartesian 
grids. An extension of our algorithm on block structured adaptive meshes is un- 
der development, and requires only modifications to the inner multigrid linear 
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solver. 

The paper is organized as follows. In Sec. 2 the governing equations for cos- 
mological radiation hydrodynamics are presented under two different assump- 
tions about the radiation-matter coupling: a two temperature model assuming 
local thermodynamic equilibrium (LTE). and a non-LTE ionization kinetics mul- 
tispecies model. In Sec. 3 we describe our solution procedures for splitting off 
the hydrodynamic calculation, and our coupled implicit radiation-ionization-gas 
energy system. Results from solution verification tests are presented in Sec. 4, 
as well as parallel scalability tests. Conclusions follow in Sec. 5. 



2. Cosmological Radiation-Hydrodynamics-Ionization Model 

We consider the coupled system of partial differential equations 

dtPb + -v& ■ Vpfa = --p b \7 ■ v b , (1) 
a a 

<9 t v b + i(v b -V)v b = --v fc -— VJ5--V0, (2) 
a a ap b a 

d t e+-v b ■ Ve = -— e — V • (pv b ) - -v b ■ V0 + G - A (3) 

a a ap b a 



1. 



d t ni + -V • (njV b ) = aijUeiij - n^rf , i = l,...,N s (4) 

d t E + • (Ev b ) = V • (DVE) - m^-E + Amr) - ckE. (5) 

Here, the first three equations ([J)-© correspond to the equations of ideal gas 
dynamics in a coordinate system that is comoving with the expanding universe 



14[. These equations correspond to mass, momentum and energy conservation, 
respectively, in which v b = a(i)x is the proper peculiar baryonic velocity, p is 
the proper pressure, and the total gas energy per unit mass is given by e. The 
modified gravitational potential is given by <f>, which satisfies the comoving form 
of Poisson's equation, 

V 2 </>= -H( Pb + Pdm - (p)), (6) 
a 

where g provides the gravitational constant, p b and pdm are the baryonic and 
dark matter densities, respectively, and (p) is the cosmic mean density. The 
densities pi are comoving, relating to the proper densities through the relation 
Pi = /0i jPrO pero(i) 3 - Here a(t) = (1 + denotes the cosmological expansion 
parameter for a smooth homogeneous background, where the redshift z is a 
function of time only. All spatial derivatives are taken with respect to the 
comoving position x = r/a(t). The hydrodynamics equations are closed as 
usual with the ideal gas equation of state, 

P 1 ,__ ,2 



P&(7-1) 2 



oi^r, (7) 
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where 7 is the ratio of specific heats, taken to be 5/3 throughout this work. 

The hydrodynamics equations are coupled with the elemental species rate 
equations (HJ) and an equation describing the flux-limited diffusion (FLD) ap- 



proximation of radiation transport in a cosmological medium ([5]) 
these equations i denotes the i th chemical species (out of N s total), n, is the 
comoving number density, n e is the electron number density, nj corresponds to 
other ions that react with the species i, and Ojj are the rate coefficients defin- 
ing these interactions [ltj [rH | ; E corresponds to the comoving radiation energy 
density. The parameter m controls whether E is monochromatic at a specified 
frequency v (m = 0), or an integrated grey radiation energy density (m = 1). 

The baryonic gas is coupled to collisionless dark matter solely through their 
self-consistent gravitational field via ([6]). The dark matter density is evolved 
using the Particle-Mesh method described in [l8, 2(J. As the N-body method 
is standard and not the focus of this paper, we do not elaborate on it here. 

2.1. Model Coupling 

In addition to the advective components of the chemistry and radiation equa- 
tions, coupling between these equations arise through a number of spatially-local 
reaction terms. The radiation energy density and chemical number densities af- 
fect the gas energy through the heating and cooling rates G and A, respectively. 
The ionization and recombination rates aj j and cmissivity 77 depend on the gas 
temperature, 

r = ( 7 -i)^r £ , (8) 

Pb Kb 

where m p corresponds to the mass of a proton, fi corresponds to the local 
molecular weight, and kb is Boltzmann's constant. Finally, the photoionization 
rate T? h depends on the radiation energy density, and the opacity k depends on 
the state of chemical ionization. 

In determining these coupling terms we distinguish between two cases: those 
in local thermodynamic equilibrium (LTE) and those that are not (nLTE). In 
the LTE case the chemical species are assumed to be in equilibrium, and hence 
their equations Q may be omitted from the time-dependent system (dJ-©. For 
problems in this regime, the coupling terms resemble those typically encountered 



in radiation- hydrodynamics simulations [12l. l2ll. \22\ 



Vlte(T) = n P B = T , 

7T 

G L TE{Pb,E) = —E, (g) 
4-7T 

K LTE {pb,T) = — tjlte(T), 
Pb 

where c is the speed of light, <jsb is the Stefan-Boltzmann constant, and Kp 
and k correspond to the problem-dependent Planck mean and total opacities 
for the gas. 

For simulations that may not be approximated as being in local thermo- 
dynamic equilibrium, these coupling terms involve the dynamically-changing 
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chemical ionization states. Here, the combined opacity depends on the local 
ionization states n,, the emissivity 77 depends on both T and n,, the gas heating 
rate G depends on both E and n.j , and the gas cooling rate A depends on T and 
nj, with the corresponding formulas given in the references 1|| 23|, 24, 25, 26[. 



2.2. Cosmological Flux- Limited Radiation Diffusion Model Details 

We derive the cosmological flux-limited radiation diffusion equation ([5]) from 
the general multi- frequency version (ljj . 

d t E v + -V • (E v v b ) = V • (DVE V ) + v-d v E v + 4^ - ck v E v . (10) 
a a 

Through assumption of a given radiation frequency spectrum, xe(v), the frequency- 
dependent radiation energy density may be written in the form E v (x.,t,v) = 
S(x, t) xe(v)- With this, we define the single "grey" radiation energy density 
used in the model H])-© as 

/>oc />oo 

E(x,t)= E v (x,t,v)du = E(x,t) XE{v)dv. (11) 

The radiation equation §5§ may then be derived through integration of the 
equation IjlOp over frequencies ranging from the ionization threshold of Hydrogen 
(hi/Q = 13.6 eV) to infinity; integration of the term v~d v E v gives rise to the 
term — - in ([SJ). We note that this approximation (JTTJ is valid only if the 
assumed spectrum xe{v) is defined such that the indefinite integral exists, as 
is the case for quasar and stellar type spectra where it scales with frequency 
as E v oc v~Pi where (3 q > 1. However, through this formulation we may also 
consider problems involving a monochromatic radiation energy density, since 
such energy densities may also be expanded as E„(pt, t, v) = £?(x, t) Xe(v) 7 where 
for radiation at the monochromatic frequency Vk, the assumed spectrum is given 
through the Dirac-delta function xe{v) — 6 Vh (y). In such cases, the term 
v-d v E v vanishes, giving rise to the parameter m in ((SJ). For standard grey 
radiation approximations, we assume a radiation spectrum of the form of either 

a power law, Xe{v) — 1 a < ~1> or as a blackbody spectrum, xe{ v ) = 

87 r/ l (|) 3 /(exp(-^)-l). 

As is standard with FLD approximations to radiation transfer, one must 
take special care in construction of the diffusion coefficient function D. In 
its simplest form, D may be written as D = 3^7, where kt = Ka + is 
the total extinction coefficient, ka corresponds to total absorption (here taken 



to be the opacity k) and ks corresponds to scattering [2l\ . Use of this form 
for the diffusion coefficient, however, results in an infinite signal speed for the 
radiative flux — S7E. To preserve causality, the analytic form of D is modified 
with a dimensionless flux-limiter whose particular form may be tuned to the 
specific problem of interest, but whose overriding purpose is to guarantee that 
the radiation transfer equation ([5]) gives the correct numerical behavior in the 
limiting cases of (nearly) isotropic and free-streaming radiation. Several choices 
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for flux- limited forms of D have been proposed in the literature [2?], [28| . We 
consider the diffusion coefficient to be of the form 



D(E) 



where 



Di(E) 










D 2 (E) 







D 3 (E) 



(12) 



c(2k t + Ri) 



6k^ 



:'./. ; /.', + Rf 



with Ri = \diE\/E, i = 1^2, 3. We note that this function has been reformulated 
from its original version [12| to allow increased numerical precision for scattering- 
free simulations involving extremely small opacities (i.e. kt = ha — K <C 1), as 
is typical in cosmology applications. 



3. Solving the Coupled System 

3.1. Operator- Split Hydrodynamics with Radiative Feedback 

Since typical astrophysical and cosmological simulations involve the hydro- 
dynamic motion of gases encountering shocks, whereas radiation diffusion and 
chemical kinetics processes are more of reaction-diffusion type, we choose to 
solve the coupled system ([T])-© in an operator-split fashion. In this approach, 
a time step t n to t n+l is taken using the general steps 

(i) Deposit the dark matter particles onto the mesh to calculate the p2 m - 

(ii) Solve for the gravitational potential <f> resulting from the densities p\, and 
Pdm using equation 

(iii) Evolve the dark matter particles using the Particle Mesh Method 0, 0, 

(iv) Evolve the hydrodynamics equations (HJ-© with a high-order, explicit- 
time upwind method. In this step, use the velocity field v^ to advect both 
the chemical number densities and radiation energy density E. 

(v) Using a high-order implicit-time method, solve a coupled reaction-diffusion 
system to obtain the time-evolved number densities n^, radiation energy 
density E and gas energy e. 

In order to allow us to split the equations CO)-© into the two steps (iv) and (v) 
above, we consider the gas energy as consisting of two components, e = e^ + e c , 
where is the fluid energy arising from the hydrodynamic evolution of the 
system, and e c is the gas energy correction arising from the couplings with 



G 



radiation and chemistry. Under this decomposition, the energy conservation 
equation ([3]) may be equivalently written as 



d t (e h + e c ) + -v b - V(e fe + e c ) = (13) 

- — (e h + e c ) - —V • (pv b ) - -v 6 ■ + G - A. 
a apf, a 

Under this splitting, the hydrodynamic solver used in step (iv) of the operator- 
split algorithm solves the system of equations 

11 

OtPb + -v 6 • Vp b = — p b V ■ v b , (14) 
a a 

<9 t v fc + - (v 6 • V) v & = --y b - — Vp - -V<f>, (15) 
a a ap b a 

d t e h + -v b ■ Ve h = -~—e h • (pv b ) - -v b ■ V0 (16) 

a a ap b a 

d t n t + -V • Kv b ) = 0, (17) 
d t E + -V ■ (Ev b ) = 0, (18) 

to evolve the solution at t n , (p% , , e™, n™, E n ), to the time-updated variables at 
t n+1 , (/Ob +1 ,Vb +1 ,e™ +1 ), and the adverted variables (n*,£*). For this step, we 
employ the Piecewise Parabolic Method (PPM) [29j , on a regular finite- volume 
spatial grid, implemented in the community astrophysics code Enzo 0, [2(], [3(| ■ 
The remainder of the coupled system, 

d t e c = ~— e c + G-A, (19) 
a 

dtrii = aijTL e nj - rnYf 1 , (20) 

d t E = V ■ (DWE) - m-E + Am] - ckE, (21) 
a 

is then solved using a fully implicit nonlinear solution approach to evolve the 
advected variables (0, n* , E* ) to the time-evolved quantities (e™ +1 , n™ +1 , E n+1 ) . 
Here we may assume that e™ = since the hydrodynamic solver uses the full 
energy at t n in its evolution, i.e. eJJ = e n . Once this step is finished, we 
compute the time-evolved total energy as the sum of the hydrodynamic portion 
e/j and the adjustments due to radiative feedback e c , i.e. e n+1 = e^ +1 + e™ +1 . 
The treatment of the implicit radiation, chemical ionization, and gas energy 
feedback system (fTT))) - ([2"Tj) serves as our focus for the remainder of this section. 
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3.2. Solving the Radiation, Ionization and Energy Feedback System 

Under a method-of-lines approach, we consider a two level, up-to-second 
order accurate theta-scheme for implicit integration of our system (Til?)) - ([2"l"1) , 

e' c l+1 + At6C n e +1 = e n c + At{6 - (22) 
n™ +1 + AteC^ 1 = < + At{6 - , (23) 
E n+1 + At9 [D n E +1 + C E +1 ] = E n + At{6 - 1) [V% + C n E ] . (24) 

Here, the parameter 6 defines the implicit integration method: 6=1 corresponds 
to a first-order implicit Eulcr method, 6=0.5 gives a second-order time-centered 
approach (i.e. Crank-Nicolson). We note that in the ensuing computational re- 
sults from section|H we have typically taken 0=0.51 to provide a nearly-second- 
order time integration while avoidin g th e "ringing" traditionally associated with 
fully time-centered approaches [3ll l32j. For the above equations, we have de- 
fined the diffusive operator 

V E = V E {E, m) ee -V ■ (DVE) , (25) 

and we have defined the local "reaction" operators as 

C e =C e {e C7 E,n q ) ee — e c - G + A (26) 
a 

£ni = £n, (n,-, e c , E) ee HiVf 1 - o,. ,n, n, (27) 

C E = C E (E,e c ,n.i) ee m-E - Atttj + ckE. (28) 
a 

The equations (|22p , (|23p and (JM]) form a coupled nonlinear system of reaction- 
diffusion equations for evolution of the fluid energy correction e c , the elemental 
number densities n^, and the radiation energy density E. Denoting the vector 
of unknowns U = (e c , n^, E) T , we first define the nonlinear residual function for 
the time step t n — > t n+1 , as 

C e (U) 

f(U) = U + Ate | C m (U) | | gl | , (29) 

V E (U) + C E {U) 

where the vectors g™ are formed using the previous time-level information from 
(|2"2"|) - (f2~4")) . In order to evolve the coupled implicit system, we solve the nonlinear 
problem f(U) = for the updated vector of unknowns U n+1 . For this nonlinear 




solve, we use a globalized Inexact Newton's Method [33l |34| , in which we apply 



an iterative process for convergence toward the solution U n+1 in the following 



Given an initial guess Uo f» U(t n+1 ), we iterate toward a solution 
U n+1 satisfying \\f{U n+1 )\\ < e < 1 (we typically choose e = 1CT 7 ): 
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1. Approximately solve the linearized Newton system, \\J(Uk)Sk + 
f{Uk)\\ < <5fc, to tolerance 5k for the correction vector Sk- Here, 
J{Uk) = -gjjfiUk), and we typically choose the tolerance as 
4 = 1CT 6 ||/(L4)II- 

2. Update the vector of unknowns as Uk+i = Uk + XkSk, where 
Afc G (Xmini 1] is the line-search parameter [3BI |36|. 

We measure convergence of the Newton iteration with the RMS norm 

11*111 



N(N S + 2) 



(30) 



where N(N S + 2) is the number of unknowns in v (N spatial cells, N s + 2 vari- 
ables), since such a norm does not grow artificially larger with mesh refinement. 
The key to efficiency of the inexact Newton algorithm lies in a fast and robust 
solver for the linear systems JS = —f. Once such a solver has been provided, 
the algorithm exhibits very fast convergence - superlinear for this choice of Sk 
33j, |37|. Moreover, for diffusive PDE systems similar to the one solved here, the 
Newton convergence rate has been shown to be independent of spatial resolution 
(3^1 , suggesting that this entire implicit algorithm should allow scalability to the 
limits of the inner linear solver. 



3.2.1. Linear Solver 

In solving the system (f2"9"f , we make one approximation within the Newton 
system matrices J(U), wherein we lag the n, dependence of Pg in (|25| to the 
previous Newton iterate. Mathematically, this results in a full Newton step for 
all but the limiter's dependence on the chemical opacities, which are instead 
converged through a fixed-point iteration. The resulting solution retains the 
accuracy and stability of the full Newton iteration, albeit with theoretically 
slower convergence. However, in practice we have not noticed any increase in 
nonlinear iterations due to this approximation, and most importantly it results 
in inexact Newton matrices with the form 



J(U) = i + Ate 



J, 



J, 
Jn 



Je,E 
Jn,E 

Je,e 



(31) 



where nearly all of the blocks are given by the spatially- local components, 

Je,e = [d e £e] Je,n = [<9 ni £ e <9 n2 C e ..] J e ,E = [t^E^e] 



<7e,, 









' d C 

^ni i -'ni 


^112 ^ni ■ ■ ■ 




' <9 B £ ni 




d C 




^ni '-'112 


^112 ^112 ' ■ ' 


Jn.E = 


9e^-H2 



[d e Cl 



JE,n = [d^C-E d n2 CE 



(32) 
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and the only block containing spatial couplings is Je,e = [9e (J>e + £e)\- Thus, 
although the Jacobian matrix contains couplings both within and between vari- 
ables, it has a very desirable structure: all inter-variable couplings occur locally 
in space, and the only nonlocal couplings arc within the block Je,e, consisting 
of a scalar-valued reaction-diffusion operator. 

In keeping with a block-structured view of the Jacobian (|3"Tj) . we rewrite the 
Newton system JS = — f in the form 



" M 


U ' 


( s M > 




( fu \ 


L 


D 




)=- 


\ fE J 



where 



M = I + Ate 



Je,e 
Jn.e 



J e,n 

L = Ate [ J E ,e J E ,n ] , 



u = Ate 



Je,E 
Jn,E 

D = I + Ate [ J e ,e 



(33) 



(34) 



SM 



[s e ,s n ] T and Jm = [fe,fn] T - We note that the only matrix containing 
spatial dependencies is D, so under an appropriate variable ordering the other 
sub-matrices are block diagonal. Hence, we may efficiently invert M to obtain 
Sm as a function of se- 

Ms M + Us E = -fM => s M = -M- 1 (f M + Us E ). 

Inserting this into the second row, we have the single equation for se, 

(D - LM- l U)s E = -f E + LM- l f M - 

Therefore, this Schur complement formulation [39| for solution of the linear 
Newton system ([3"3")) proceeds with the following steps: 

(i) Set f u = M~ l f M . 

(ii) Solve the system (D ~ LA1~ 1 U)se = —f E + L/m for s E . 

(hi) Recover the remaining solution pieces, sm = — /m — M~ 1 Us E . 

We examine each of these steps below. 

The step (i) corresponds to solving the linear system 



/ + Ate J e ,e 

Atej n . e 



At6»J e , n 
/ + Ate J n , n 



fe 
fn 



fe 
fn 



(35) 



Due to the spatial locality of each component in M, we order the equations and 
unknowns in this system so that application of M, and more notably M _1 , may 
be performed independently in every spatial cell. Such solves consist of dense 
matrix algebra on (l + iV s ) x (l + iV s ) linear systems (for N s chemical densities). 
In addition, at this step we compute the matrix M^ 1 U through one additional 
solve with M, which will be used in the following steps. 
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The step (ii) corresponds to solving the system {D—LAI~ 1t I)se = Lfu—fB- 
This is denoted the Schur complement system, with the matrix S = D—LM~ l U. 
We note that due to the spatially-local nature of the matrices L and M~ 1 U, 
we may form S by constructing the diffusive sub-matrix £>, followed by updates 
to the diagonal entries corresponding to the entries of LM~ 1 U. Similarly con- 
struction of the right-hand side L/m — Je may occur independently at each 
spatial location. Once this system has been computed, we use a multigrid- 
preconditioned conjugate gradient parallel linear solver from the HYP RE library 



[40|, |41| to perform the scalar-valued solve, se = 5 _1 (i/A/ — /e). We note that 
this is the only step in the solution of the Jacobian systems ([55)1 that requires 
communication between processors. Moreover, we point out that in recent tests 
the HYP RE library has demonstrated ideal weak scaling up to over 100,000 pro- 
cessors for diffusion problems similar to the one encountered in this work [42[. 
As this solver comprises the majority of the non-local components within our 
nonlinear solver, we therefore expect similar scalability for the overall implicit 
solution approach described here. 

The final step (iii) in solution of the system (|33[) is to recover the solution 
components % = (s e , s n ) T via the system sm = — fit — M~ 1 Use- Again, since 
we have already computed the spatially- local matrix M~ lT J and the vector f^j = 
M Jm m step (i), we may trivially obtain the remaining solution components 
through cell-local matrix- vector products and vector operations, s M = —/m — 
(M- l U)s E - 

3.2.2. Multiphysics/ 'Cosmology Units 

As with any multi-physics system, special care must be taken when solving 
such systems computationally due to disparate scales between variables and 
equations. This problem is especially evident in cosmology applications, where 
in CGS units one may typically enounter specific gas energies on the order of 
10 12 , number densities on the order of 10 -27 , and radiation energy densities 
on the order of 10 -15 , with all proper density values decreasing in time due to 
cosmological expansion. To this end, we define the scaled variables 

e c = e c /u ei E = E/u E , a i =n i /u n , (36) 
x = x/u x , t = t/ut, 

where the constants u e , ue, u n , u x and u t correspond to the typical magnitudes 
of gas energy, radiation energy density, chemical number density, length and 
time at the start of the simulation. We note that due to our use of comoving 
values for E, and x, these constants are all redshift-independent, with the 
proper values of these quantities given by 

^proper = E / a 3 (t) = E 



x a(t) = x u x a(t). 



Improper = ^i/o?{t) = hi-^—, (37) 

a (t) 
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The constants are supplied on a problem-dependent basis, to allow for adapt- 
able, on-the-fly rescaling of simulations ranging from normalized test problems 
to cosmological reionization. With these rescaled variables, we rewrite our equa- 
tions (|2"2"f - (f2~Xf as the normalized system 



gn+l 



At8£™ +1 



At(9 - l)C n e , 



E 



Ate v n E +1 



rn+l 
E 



E n + At{6 - 1) 



R 



(38) 
(39) 
(40) 



Here the operators £ e , C n , Ce and T>e have correspondingly absorbed the renor- 
malization constants u*. These equations, along with the normalized solution 
vector U = (e C) ni,E) T are then used within the solution strategy described in 
section [~ 



3.2.3. Adaptive Time Step Selection 

A strong appeal of using implicit methods is their stability with respect to 
time step size; however such freedom gives rise to the qucstionof what time step 
should be used. At one extreme, we may choose a large step to achieve overall 
efficiency of the simulation, with little to no knowledge of the resulting temporal 
accuracy. At the other extreme, we may choose a very small time step for an 
accurate solution, resulting in inefficient simulations due to the increased cost 
of solving the nonlinear systems at each step. As the approach described here is 
operator-split, in which the hydrodynamics is solved using an explicit approach, 
we are therefore bound by the hydrodynamic CFL stability limit; however for 
most problems involving radiation and chemical ionization, the dynamic time 
scales of interest remain significantly faster than the hydrodynamic time scale. 
Thus the question of how to adaptively choose the time step size remains. 

To that end, we adaptively choose the time steps as the largest possible 
that additionally satisfy a prescribed accuracy requirement. We estimate this 
accuracy through comparison of the updated solution U n+l with an explicit 
predictor for that solution U pred . Defining the weighting vector in a spatial cell 
i for the variable v as 

(41) 



Trn+lj-Tpredi 



1. 



(which assumes normalized values oiU v ), we estimate the local accuracy of the 
current time step as 



1 



N(N S + 2) 



jjn+l jjpred 



1/P 



(42) 



where we have used the standard p-norm (including p = oo as the 'max' norm, 
in which case we do not divide by N(N S + 2)), and where the quotient inside 
the norm is taken pointwise. With this estimate, we set the new time step to 



At n+1 



Ttol Af 

£loc 



(43) 



12 



which should provide the maximal value that still satisfies the desired integration 
accuracy tolerance, T to i, assuming that U pred approximates the time-evolved 
solution U n+1 to O(At). 

Here, the vector cu is designed so that ei oc estimates the average relative 
change in each solution component, and includes the harmonic mean of the pre- 
dicted and new states to allow increased robustness in the case of cosmology-type 
problems where variables change by orders of magnitude across cells and time 
steps. The value of p is typically taken to be oo in the ensuing test problems; 
however such a choice may limit parallel scaling since such a measure is sensi- 
tive to pointwise changes, of which there are many more as dynamics propagate 
throughout an increasingly refined domain. Lastly, we use the explicit predic- 
tor as the initial guess for the Newton method, Uo, which we describe in the 
following section [3.2.41 

3.2.4- Explicit Predictor 

A well known property of Newton's method is that its robustness and effi- 
ciency benefit greatly from an accurate initial guess. To this end, we provide 
the predicted initial Newton iterate 



i.e. we use an initial guess given by the 0(A£)-accurate explicit Euler update 
to the coupled system ([T^|) - ([2T|) . As this provides only an initial guess to the 
solution, its instability at larger At will not affect the temporal stability of the 
overall method, since the solution to each step must satisfy the implicit system 
(|2"2"|) - (f24"]) . However, as we use an adaptive time-stepping strategy, for very fast 
dynamics (that give rise to very small At), such an initial guess may already 
satisfy the nonlinear tolerance ||/(J7o)|| < £ and the solver will not require any 
Newton iterations, effectively allowing an adaptive explicit/implicit simulation 
of the coupled system (Til)]) - ([2"Tj) . 

3.2.5. Adaptive Computation with Supplied Radiation Spectrum Xe{v) 

The final detail that we describe in this solution approach relates to the 
choice of assumed radiation spectrum xe(v). As noted in section \2.2\ we may 
choose cither a monochromatic or an integrated "grey" radiation equation, based 
on the choice of this assumed spectrum. This choice affects all terms involving 
the radiation energy density E v in the general radiation energy equation pop 
and in the coupling terms G and r? . As each of these terms involve a product 
of the form f{v)E v , integration over v converts these to 



We therefore allow a user-defined functional form for xe{v), which we then 
numerically integrate to high accuracy upon initialization of the simulation, 




(44) 




(45) 
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providing the relevant constants necessary to convert the ^-dependent equation 
(fT0|) to the monochromatic or grey integrated equation §5§. 

4. Numerical Results 

Wc present test problems designed to verify the accuracy of the radiation 
diffusion and chemical ionization modules in conjunction with hydrodynamical 
fluid motions. Since a number of distinct processes may compete for impor- 
tance in a full simulation, we begin with simple tests that isolate various single 
components, and subsequently build upon those results with more sophisticated 
problems that couple additional physics. We begin with a radiation diffusion 
problem ( ^4.1|) that exercises the diffusion term of the radiation equation ||5j) 
in the absence of energy, chemical, or hydrodynamic coupling; this test is fol- 
lowed by an examination of the matter-radiation coupling terms ( H4.2|) in an 
infinite uniform medium, a diffusion wave with material coupling ( ^4.3[) . and a 
non-equilibrium radiating shock problem f H4.4[) . all of which assume chemical 
equilibrium. Our attention then turns to problems including ionization, be- 
ginning with an Hll ionization front propagating through a static, isothermal 
medium ( H4.5I) . followed by a cosmological 1- front propagation problem that ex- 
ercises the cosmology terms and units f ^4.6|) . We then consider a fully-coupled 
radiation- hydrodynamics- ionization calculation 7[) . This section concludes 
with additional calculations f N4.8p demonstrating the parallel scalability of the 
radiation diffusion module, which of all components places the highest commu- 
nication demands upon a domain-decomposed parallel calculation. 

We note that for all test problems except ( §4.6|) and ( §4.8j) . we use a non- 
cosmological problem (i.e. z = and a = 1). In problems ( §4.6j) and ( §4.8|) . the 
cosmological parameters are described therein. 

4-1- Free-Streaming Radiation 

Because the standard diffusion equation is parabolic, the associated signal 
speed of the diffusion variable is formally infinite. However in reality radiation 
fronts propagate at speeds bounded by the speed of light in vacuum, so we 
modify the diffusion coefficient in our radiation energy equation ([3]) with a flux- 
limiter, as discussed in section 12.21 Our first test problem verifies the correct 
action of this limiter by examining the propagation of a planar radiation front 
through a transparent medium. Radiation is assumed to propagate along the 
x-axis of our computational mesh; a Dirichlct boundary condition is imposed 
on the left boundary specifying an incident radiation energy density of 1.0 erg 
cm~ 3 . Physically, the expectation is that with a sufficiently small (but nonzero, 
due to numerical constraints) opacity, a sharp radiation front will move through 
the domain at the speed of light. The Planck and Rosseland mean opacities 
are assigned a constant value of 10~ 6 cm -1 , ensuring an essentially transparent 
medium. The spatially uniform initial value of the radiation energy density is 
10~ 4 erg cm~ 3 . 

The computational mesh has a domain length of 1.0 cm along the propaga- 
tion direction of the light wave. We have run the problem for 8.3391 picoseconds, 
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Figure 1: Curves of E vs. x for mesh sizes of 128 (red), 256 (orange), 512 (blue), 1024 (green), 
and 2048 (violet) zones. The analytical solution (black dashed line) is a step function centered 
at x = 0.25cm. 

which is one quarter of the light-crossing time for this length. Figure [T] shows a 
series of curves resulting from calculations at mesh sizes of 128, 256, 512, 1024, 
and 2048 zones along the x-axis. The dashed line indicates the expected loca- 
tion of the radiation front, ct, where c is the speed of light and t the evolution 
time. In the absence of the flux limiter, the numerical curves would give the 
formal t— >oo solution for the diffusion equation, which for our problem param- 
eters would be a nearly horizontal profile for E throughout the domain (given 
the nearly zero opacity). That the curves capture the correct location of the 
radiation front is due entirely to the action of the limiter. The sharpening of 
this front with increased resolution is evident. 

Also apparent is a slight lag (about 0.01 cm) between the analytical loca- 
tion of the radiation front and the numerical location, taken as the common 
intersection point of the numerical curves. The size of this lag depends upon 
the choice of the adaptive time step. In the language of f ^3.2.3|) . we compute 
At using p = oo and we vary r to i, which here corresponds to the maximum 
allowed fractional change in the radiation energy density per timestep. Figure[2] 
illustrates this effect by showing 3 curve pairs, each of which has been com- 
puted with a different choice of T to i: 0.01 (green), 0.05 (blue), and 0.25 (red). In 
each colored pair, the solid curve was obtained using 512 zones, and the dashed 
curve shows the 128-zone result. Two curves at different mesh resolutions are 
provided to identify the "consensus" value of the light front location via their 
point of intersection. This location is seen to converge as r to i — > 0. 
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Figure 2: Radiation energy profiles at 128 (dashed curves) and 512 (solid curves) zones for r to i 
of 0.01 (green), 0.05 (blue), and 0.25 (red). Intersections of associated 128 and 512-zone curves 
indicate the "consensus" location of the radiation front for a given choice of r to i . Convergence 
to the correct value is readily observed. 

4-2. Matter- Radiation Equilibration in a Homogeneous Medium 

We now consider a problem in which e and E are spatially uniform but 
are initialized to values far away from equilibrium. This case thus isolates the 
matter-radiation coupling terms in the gas and radiation energy equations. The 
parameters for this test were published by Turner and Stone [43( , who assumed 
an isotropic medium characterized by a single opacity of 4 x 1CP 8 cm -1 , a 
gas density of 10" 7 g cm -3 , and an average particle mass of 0.6 m#, where 
vciH is the mass of a hydrogen atom. Coupled to this medium is a radiation 
field with a uniform value of 10 12 erg cm" 3 . From this value we compute a 
"radiation temperature", T r = (E/ar) 1 / 4 , of about 3.4 x 10 6 K. Here we have 
defined a r as the radiation constant, 7.56 x 10" 15 erg cm" 3 K -1 . Two cases 
are considered: one in which the initial gas energy density is e = 10 10 erg 
cm" 3 , and one in which the initial value is e = 10 2 erg cm -3 . For the stated 
parameters, these energies correspond to gas temperatures of roughly 4.8 x 10 8 
and 4.8 K, respectively, which therefore bracket the radiation temperature. In 
both cases, however, the initial radiation temperature is sufficiently high that 
the radiation energy density should remain constant to good approximation as 
the gas evolves to thermal equilibrium. To see this clearly, consider the effective 
heat capacity of a unit volume of the radiation "gas" as compared to that for the 
material. For radiation, this number is simply 1.0 cm x a r T 3 , which evaluates 
to roughly 3 x 10 5 erg K . In contrast, the gamma value and mean particle 
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Figure 3: Evolution to thermal equilibrium of a medium with an initially high (solid line) 
and low (dashed line) gas energy density. The analytical equilibrium value is shown by the 
horizontal dotted line. 



mass translate to a specific heat (erg g _1 K _1 ) of roughly 2.0 x 10 8 , which yields, 
for our assumed density, a material heat capacity of 20 erg K _1 for a unit control 
volume. The physical result is that the radiation field has an effectively infinite 
thermal reservoir when compared to the material. 

If the radiation energy density is formally assumed to be constant, the gas 
energy equation may be written as a simple ODE: 

e = ckE - 4irKB(e), (46) 

where B is the temperature-dependent Planck function 

b{t)= c^ t4 m 

Using the ideal-gas law ([7|) , we write B as a function of e and solve the simplified 
gas energy equation for the equilibrium value of e such that e = 0, 



3 / P k B \ ( 1/4 



-'eq 



2 \0.6m H 



(48) 



Notice that this expression is nothing more than the ideal gas formula for e 
evaluated at the fixed radiation temperature, the expected result. 

The results of our two test calculations are shown in figure [3] Both tests 
were run in a small box domain (4 3 zones) with triply-periodic boundaries. The 
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io- v 
10- 10 
io- 4 
io- 7 
io- 7 


io- lJ 
io- 9 
io- 9 
io- 12 

IO -3 


7.04 x 10~ 12 
3.07 x 10~ 13 
3.55 x 10 -11 
7.04 x IO -12 
7.04 x IO -12 


1.37 x IO -12 
3.45 x IO -14 
1.83 x IO -11 
1.37 x IO -12 
1.37 x IO -12 



Table 1: Conservation of total energy for the matter-radiation equilibration test. Relative error 
in energy conservation H49H for both low-to-high and high-to-low temperature equilibration, 
for various nonlinear and linear solver tolerances e and 8, respectively. 

case of T(0) > T r is indicated by the solid curve; the low-T case is shown by 
the dashed line. The horizontal dotted line has been placed at e = e cq . Both 
energy curves converge to the correct result. Note that in this test the opacity 
serves only to control the timescale to reach thermal equilibrium; neither the 
value of the equilibrium energy nor the validity of our assumption of constant 
E are dependent on the value of k. While figure [3] demonstrates convergence 
to the correct asymptotic value of the gas energy, it provides no information as 
to whether the rate at which it approaches this value is correct. A quantitative 
assessment of this latter metric is provided by our next test problem. 

Wc further use this test to examine the conservation properties of the coupled 
radiation and gas energy solver. In Table [T] we show the value of 

J \E total (t) - E total (0)\ dx 

jE total (0)dx 1 ) 

for both tests at the final time i=2.5e-7 sec, run using a variety of nonlinear 
and linear solver tolerances, e and 5, respectively. We note that in all cases, the 
total energy is conserved to more than 10 digits of accuracy. Moreover, while 
the conservation is weakly dependent on the nonlinear solver tolerance, it is 
entirely independent of the linear solver tolerance. This behavior is most likely 
due to use of the Schur complement formulation f £|3.2.ip . that exactly solves 
for coupling between variables to floating point roundoff, leaving the iterative 
linear solver to handle only the radiation equation. Wc further note that this 
is an ideal problem to test conservation of the coupled solver, since it is the 
only test considered that uses a closed system. We further comment that since 
the PPM finite-volume method is constructed to satisfy conservation, and the 
implicit subsolver achieves conservation to high accuracy, overall conservation 
of the coupled solver follows. However, we note that for problems utilizing non- 
periodic boundary conditions, chemical ionization cooling, gravitational heating, 
or cosmological expansion, the model no longer represents a closed system and 
therefore will not conserve energy. 

4-3. Non- Equilibrium Marshak Waves 

This test exercises both radiation diffusion and the physics of matter-radiation 
coupling. Non-equilibrium Marshak waves characterize the evolution of the ra- 
diation field in an initially cold, uniform halfspace on which a radiation source 
is imposed. The particular form of the Marshak problem described here is 
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originally due to Pomraning 44|. The problem was re-examined by Su and 
Olson [451 ] . who derived semi-analytic exact solutions for the radiation and gas 
energy densities and tabulated select values of them on a grid of space and time 
values. The problem considers a formally 1-D semi-infinite domain in which 
(z, t) denote dimensional space and time coordinates. Ignoring hydrodynamic 
motions, Su and Olson write simplified forms for the radiation and material 
energy equations as 

d t E{z, t) - d z [^d z E(z, i)] = ck [a r T 4 (z, t) - E(z, t)] , (50) 
c v (T)d t T(z,t) =c K [E(z,t)~a r T 4 (z,t)] , (51) 

in which k is the constant opacity, a r the radiation constant as defined previ- 
ously, and c v the specific heat of the material. Note that the flux-divergence 
term in (|50[> assumes pure diffusion with no flux limiter. The matter tempera- 
ture, T, is assumed to be related to the gas energy, e, via 

e = pc v T. (52) 

As written, ([50)) and (fSTj) are coupled nonlinear PDEs in the dependent variables 
E and T. Su and Olson linearized the equations by choosing the following form 
for the specific heat: 

c v = aT 3 , (53) 

where a is an arbitrary constant. The T 3 dependence of c„ on temperature has 
two effects: it allows equations ([50]) and (jSTj) to be written as linear ODEs in 
E and T 4 , and it gives the heat capacity of the material the same temperature 
dependence as the effective heat capacity of the radiation field. With an appro- 
priate choice of a, a problem can therefore be designed in which the material 
and radiation will both evolve significantly in space and time. 

The description of the problem is completed with a specification of the 
boundary conditions. A Marshak boundary condition is applied to E at z = 0: 

E(0,t) - ^-8 z E(0,t) = -F inc , (54) 

OK C 

where -Fine is the incident flux at z = 0. The boundary condition at z = oo, and 
initial conditions at t = are 

E(oo,t) = 0, E(z,0) = T(z,0) = 0. (55) 

Su and Olson construct the linearized equations by defining dimensionless 
independent and dependent variables, (X, t) and (u, v) such that 

X = zkV3 7 u(X,t) = (^-^jE(z,t), (56) 

'4et r CK N 



a J v ' V4Fi 



v(X,r) = \-^\T\z,t) 
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Figure 4: Curves of e (dashed line) and E (solid line) vs. dimensionlcss distance, at dimen- 
sionless times r = 1, 3, 10, 30, and 100 (curves shift upward as t increases). Analytic values 
for e (squares) and E (circles) are shown for r = 1, 3, and 10; values for E are also provided 
at t — 100, by which time e and E are nearly equilibrated. 



With these definitions, and letting e = 4a r /a, equations (f50|) . (f5"Tj) . and ([54)) - 
(|B"5")) become 

ed T u{X,r) - 8% 2 u(X,t) = v(X,t) - u(X,t), (57) 
d T v(X,r) = u{X,t) - v(X,t), (58) 

tt(0,T) - -^flbr«(0,r) = 1, (59) 

u(oo,r) = 0, (60) 
u(X,0) = v(X,0) = 0. (61) 

The Marshak boundary condition represented by (|59|) enforces the constraint 
of constant flux on the left boundary. This is an example of a "mixed" or Robin 
boundary condition, and as such requires special treatment in Enzo. For the 
purposes of this verification test, we implement this boundary condition by 
imposing a Dirichlet condition with a time-varying value of u computed from 
[45|'s equation 36, evaluated at X = 0. Because the integrands in their equation 
are highly oscillatory for r« 1, we substitute the asymptotic expression given 
by their equation 51 when r < 10~ 5 . 

Figure [4] shows results from a high-resolution simulation with 2048 zones 
along the X coordinate. The exact solution values tabulated by [45| span the 
range < X < 10. Since the right boundary condition is specified at X = oo, 
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Figure 5: Marshak problem convergence: curves of u (solid lines) and v (dashed lines) at 
mesh resolutions of 128 (red), 256 (orange), 512 (blue), 1024 (green) and 2048 (black) zones. 
Reference solution values are indicated by open circles and squares; the halving of relative 
error with each doubling of mesh size is readily apparent. 

we choose our domain X <G [0, L] such that L is sufficiently large (about 35) 
for the evolution time of interest that the Dirichlet condition X(L) = may be 
reasonably applied. We choose opacity and coupling parameters k = e = 1.0 
cm -1 . The curves indicate profiles of u (dashed lines) and v (solid lines) for 
r values of 1, 3, 10, 30, 100. The squares and circles indicate exact values of 
u and v, respectively. We have indicated these values on corresponding curves 
at evolution times sufficiently early that the material and radiation have not 
yet had time to equilibrate. Figure O shows a resolution study for the curves 
computed at r = 1. Curves at mesh sizes of 128 (red), 256 (orange), 512 (blue), 
1024 (green), and 2048 (black) are shown. Each calculation is performed with 
the timestep restriction T to i = 0.05. Because this treatment allows for adaptive 
timesteps, the evident first-order rate of convergence measures the combined 
effect of time and space discretization methods. 

4- 4- Subcritical Radiating Shock Waves 

Wc now add hydrodynamic motions to our mix of physics by examining the 
propagation of shock waves for which the radiation energy plays a significant 
role in the shock structure and evolution. Radiating shock waves represent a 
broad class of phenomena figuring prominently in both astrophysical and terres- 
trial applications. The particular formulation of the problem we present is due 
to Lowric and Edwards who considered the propagation of planar, steady 
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Figure 6: Subcritical radiating shock test using 4096 spatial zones. Gas and radiation tem- 
peratures are plotted in units of their preshock values: T is the solid curve; T r is the dashed 
curve. Semi-analytic values for T and T r are indicated by circles and squares, respectively. 



shock waves in the grey nonequilibrium diffusion limit. Under the assumption 
of steady flow, (4(| transform the coupled gas and radiation energy equations 
into a set of nonlinear ODEs in dimensionless gas and radiation temperature 
variables, which must be integrated numerically to achieve semi-analytic solu- 
tions. Nonetheless, their radiation diffusion model corresponds identically to 
that implemented in Enzo in the grey LTE limit, and the unique structure of 
the post-shock material temperature profile for a given Mach number makes 
this problem an excellent verification test for computer codes. 

We have run the Mach-2 test case described in (46[. The computational 
domain has a length of 0.1 cm. The material has a uniform initial density of 1.0 
g cm~ 3 , a constant specific heat of 2.218056 x 10 12 erg g _1 eV _1 , and a uniform 
initial velocity of 1.9475 x 10 5 cm s _1 . The material and radiation are assumed 
to be in thermal equilibrium at t = at a temperature of 121.6 eV. Outflow and 
reflecting boundary conditions are imposed upon the left and right boundaries, 
respectively, resulting in a shock wave that forms near the right boundary and 
propagates to the left. The total evolution time is 1.73325 nanoseconds. 

Figure [5] shows the result of a high-resolution simulation (4096 zones along 
the propagation axis). The curves represent the dimensionless gas (solid curve) 
and radiation (dashed curve) temperatures, T and Tf. The circles and squares 
are taken from exact solution data kindly provided by R. Lowrie for this param- 
eter set. Both the gas and radiation have dimensionless far-field temperatures 
of 1.0 in the pre-shock state. Examining the gas temperature curve, there are 
three features of particular interest: the precursor, in which the material is pre- 
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Figure 7: Subcritical radiating shock convergence: gas temperature vs distance from the shock 
for mesh resolutions of 128 (red), 256 (orange), 512 (blue), 1024 (green), 2048 (cyan) and 4096 
zones (black). Semi-analytic reference values are indicated by open circles. 



heated ahead of the shock front by the radiation wave which travels ahead of 
the shock; the Zel'dovich spike, shown by the overshoot in temperature at the 
shock front, and the radiation relaxation region, delineated by the decline in the 
material temperature to its eventual far-field postshock value. Letting T p de- 
note the maximum prcshock value of the gas temperature in the precursor, and 
T\ the asymptotic postshock value, we note that the property T p < T\ identifies 
this calculation as an example of a subcritical radiating shock. In the limit of 
high Mach number, T p can become equal to (but never exceed) T\, such a shock 
wave is referred to as supercritical. 

As vividly demonstrated by [4(|, the strength of the precursor, the height 
of the Zel'dovich spike, and the precise temperature structure in the relaxtion 
region are extremely sensitive to the Mach number. While the case we have 
shown is subcritical, it lies near the limit for which a multidimensional code can 
reasonably capture this structure without resorting to adaptive mesh refinement. 
The degree to which we resolve this structure as a function of resolution is 
shown in figure [3 in which we magnify the region near the shock and show 
gas temperature curves for mesh sizes of 128, 256, 512, 1024, 2048 and 4096 
zones. As shown in [Zfjj], raising the Mach number results in a dramatic increase 
in the height of the spike and narrowing of the relaxation region; a proper 
representation of the postshock structure in a supercritical shock with Enzo 
must await the implementation of adaptive mesh refinement in our radiation 
module. 

Since this problem considers coupled radiation and hydrodynamics, we also 
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examine how the adaptive time step selection strategy from Section 13.2.31 com- 
pares with the hydrodynamic CFL-limited time step. For this problem, the 
average radiation time step ranged from l.lc-5 down to 9.4e-5 for the coarsest 
(128-cell) to finest (4096-ccll) grids, exhibiting a near-constant time step selec- 
tion that tracks evolution of the radiation field. For these same problems, the 
hydrodynamic CFL limits on the time steps were 3.3e-4, 1.7e-4, 4.2e-5, 2.1e-5, 
1.0e-5 and 5.0e-6. Hence, for most problems the stiff radiation time scale lim- 
its the overall time step size, until very fine grids where the mesh-dependent 
hydrodynamic CFL stability condition becomes more restrictive. 

4-5. Isothermal Ionization of a Static Neutral Hydrogen Region 

Our first test problem incorporating ionization chemistry is due to Iliev et 
al. This problem combines radiative transfer and hydrogen ionization in a 
static astrophysical region. The physical situation of interest is the expansion 
of an ionized hydrogen (HII) region in a uniform gas around a single monochro- 
matic ionizing source emitting iV 7 = 5 x 10 48 photons per second at the ioniza- 
tion frequency of hydrogen (hv = 13.6 eV). We enforce a fixed gas temperature 
of T = 10 4 K, and a static hydrodynamic state (i.e. p\ = v& = e = 0). In such a 
problem, the radiation source should rapidly ionize the surrounding hydrogen, 
and then should develop a spherically-propagating ionization front (I-front) that 
propages quickly at first, slows, and then eventually stagnates at an equilibrium 
position referred to as the Stromgren radius, where ionization (HI— s-HII) and 
recombinations (HII— s-HI) balance. For this scenario, the analytically-provided 
I-front radius is given by 



and the recombination time is given by t rec = (a^n^) -1 . Here, as = 2.59 x 
10 -13 cm 3 s _1 is the case B hydrogen recombination coefficient. 

We have the following problem parameters: the domain size is L = 6.6 kpc 
in each direction; the initial gas number density is n# = 10 -3 cm -3 ; the initial 
radiation energy density is E = 10 -20 erg cm -3 ; the initial ionization fraction 
(HII/H) is 0.0012; the ionization source is located in the lower corner of the box 
(the (1,1,1) cell); we use reflecting boundary conditions at the x-, y- and z- 
left boundaries, and outflow conditions at the corresponding right boundaries. 
For these parameters the Stromgren radius r$ = 5.4 kpc, the recombination 
time t rec « 3.86el5 s (« 122.4 Myr), and the total simulation time is 500 Myr 
(« At rec ). The implicit solver parameters used were a convergence norm of 
p = 2, desired solution tolerance of r to i = 0.01, time-step parameter of 9 = 0.51, 
and nonlinear solver tolerance of e = 10~ 7 . 

In Figure [8] we plot the spherically-averaged I-front position and radius with 
respect to time, for various spatial mesh sizes. The I-front position is computed 
from our results as the distance at which cells transition from below 50% to 



rs[l - exp(-t/t rec )} 1/3 , where 



(62) 




(63) 
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Figure 8: Left: evolution of the I-front position, analytical solution, and relative error for 
varying resolutions. Right: computed I-front velocity and relative error. 

above 50% HII fractional density. Assuming a spherical HII region, we compute 

1/3 

this radius as r$ — (S^-) , where V is the volume comprised of all ionized 
cells (i.e. where ny/j/n^ > 0.5), and the additional factor of 8 arises due to 
the fact that our source is in the corner, so we must mirror V into the other 
7 octants. We also plot the error in the computed I-front radius and velocity 
for varying mesh sizes. As can be clearly seen, the computed I-front position is 
highly accurate, even for coarse spatial grids, with the corresponding accuracy 
increasing as the mesh is resolved. 

In Figure [9] we show cross-sections of the radiation energy density through 
the ionization source for 16 3 and 128 3 grids. We note that although the spher- 
ical front is jagged for coarse grids, as the mesh is refined we approach the 
physically-accurate spherical profile. Moreover, this demonstrates that although 
the flux limitcr Q12[) is based on one-dimensional derivatives, it does not result 
in anisotropic propagation biased along axial directions. 
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log radiation density, t = 1 00 Myr log radiation density, t = 1 00 Myr 




Figure 9: Cross-sections of the radiation density through the ionization source at t=100 Myr: 
16 3 grid (left) and 128 3 grid (right). Note the convergence to the spherical analytical shape. 



4-6. Cosmological Radiative Ionization 

We now perform the same test as above, but within a cosmofogically-expanding 
universe. The problem is originally due to Shapiro & Giroux [4jj, and combines 
cosmology, radiative transfer and chemical ionization. Here, the physics of in- 
terest is again the expansion of a HII region in uniform gas around a single 
monochromatic (hv = 13.6 cV) ionizing source. Again, the ionization front 
should propagate quickly at first, approaching the Stromgren radius, but then 
should begin to lag behind as cosmological expansion drives the Stromgren ra- 
dius outward faster than the I-front can propagate. Due to the cosmological 
expansion, the Stromgren radius changes in time, and is given by 



rs(t) 



3iV 7 



47ra_Bn ff (t) 5 



1/3 



(64) 



where although the Hydrogen number density n# is spatially static, it diminishes 
due to cosmological expansion by a factor of a~ 3 (t). Defining the parameter 
A = aB^H,%l 'Hq/ '(1 + Zi), where the subscript i refers to the quantity at the 
initial redshift zq, the analytical solution is given by 



ri(t) = r Sli [ Ae- T W / ' ' e T ^ [1 - 2q + 2q (l + *)/5] 




(65) 



where 

r(a) = A [6 q 2 (1 + z t f] ^ [F(a) - F(l)} , (66) 
F(a) = [2 - Aq - 2<7o(l + Zi)/a] [1 - 2q + 2q (l + z t )/a} 1/2 . (67) 
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Here, go is the cosmological deceleration parameter and z» is the initial redshift. 
We perform four of the tests provided in the original paper [47j |: go of 0.5 and 
0.05, and z% of 4 and 10. These correspond to the parameters: 



(q ,Zi) = (0.5,4): U w 80 kpc, H = 0.5, TO = 1.0, fi A = 0, Q b = 0.2 
(9b, Zi) = (0.05,4): L, a 60 kpc, H = 1, Sl m = 0.1, fi A = 0, tl b = 0.1 
(<?o, ^i) = (0.5, 10): Li » 36 kpc, ff = 0.5, fi m = 1.0, fi A = 0, fl b = 0.2 
(g ,*i) = (0.05, 10): L t » 27 kpc, #„ = 1, « m = 0.1, tt A = 0, tl b = 0.1 



where Li is the initial box size, Hq is the Hubble constant, £l m is the contribution 
of all non-relativistic matter to the gas energy density at z = 0, in units of the 
value required to close the universe, similarly Q\ and Q b are the contributions 
of the cosmological constant and the baryonic matter to the energy density, 
respectively. These two types of cosmology result in slightly different functions 
for the expansion coefficient a. For the case of go = 0.05, this value comes from 



cquqations (13-3) and (13-10) in [48j]. For the case go = 0.5, we use the standard 
formula a = (1 + z)^ 1 . We begin all problems with an initial radiation energy 
density of E = 10~ 35 erg cm -3 and an initial ionization fraction (HII/H) of 0. 
The initial density is dependent on go, with p b ,i = 1.175 x 10 -28 g cm -3 for 
go = 0.5, and p b .i = 2.35 x 10 - 28 g cm~ 3 for go = 0.05. All simulations are run 
from the initial redshift Zi to z = 0. All other problem parameters are identical 
to those in £14.51 All implicit solver parameters are also identical to those in 
§4.51 but with desired solution accuracy Ttoi = 0.001 and inexactness parameter 
5 k = 10- 13 ||/(C/ fc )||. 

In Figure [10] we plot the scaled, spherically-averaged I-front position with 
respect to scaled redshift for each of the four tests (with axes identical to |47j |. 
Figure la), as well as the corresponding plots for just the Zi = 4 tests along 
with their analytical solutions. These solutions all used a uniform 128 3 spatial 
mesh. In Figure [TT] we plot the error in the computed I-front radius for varying 
mesh sizes for the two cases of go = 0.5 and go = 0.05 with Zi = 4. Again, the 
accuracy in the computed I-front position improves as the mesh is resolved. 

4-7. Hydrodynamic Radiative Ionization 

We now incorporate hydrodynamic motion into the mixture of physical pro- 
cesses, and examine a problem due to Whalen & Norman that combines 
radiation, hydrodynamics and chemical ionization (but not cosmology). The 
problem is nearly identical to that from £14.51 but now in a dynamic medium 
(varying temperature, density and velocity). Again, the physics of interest is 
the expansion of a HII region in an initially uniform gas around a single ionizing 
source, though now the source emits N 7 = 5 x 10 48 photons per second with a 
frequency profile given by a Tg = 10 5 blackbody spectrum. Here, the ionization 
front should propagate quickly at first, slowing until it reaches the Stromgrcn 
radius (|63p , at which point the I-front transitions from radiation-driven (R-type) 
to dynamically-driven (D-type) , and the high pressure of the ionized and heated 
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r.(t)/r.(t) vs redshift 
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r.(t)/r (t) vs redshift, z =4 
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Figure 10: Left: I-front radii vs. redshift for the four tests, qO = 0.5 and 0.05, and z\ = 4 and 
10. Right: I-front radii vs redshift for the Z{ = 4 tests; analytical solution values are given by 
the open squares. 



Error in r(t)/r (t) vs redshift, qO - 0.5, zO = 4 
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Figure 11: Convergence of I-front radius vs. redshift for the two cases go = {0.5,0.05} and 
Zi = 4 as the mesh is refined: spatial meshes shown are 16 3 (black dotted), 32 3 (blue dashed), 
64 3 (red solid), and 128 3 (green dot-dashed). 
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gas inside the HII region continues to push the I-front out past the Stromgren 
sphere. The expansion will finally stall when a pressure equilibrium has been 
reached, at a radius Tj = (2Tj/T e ) 2 ^ 3 rs, where Tj is the temperature of the 
ionized gas behind the front, and T e is the temperature of the ionized gas ahead 
of the front. Analytical models for the initial radiation-only phase predict 



rs 



1 



-t/t r 



1/3 

(68) 



where t rec = [a^ (Tariff] 1 is the recombination time (assumed constant in this 
phase). Analytical models for the subsequent pressure-only phase predict 

7c s ^ 4/7 



. (69) 



where c s = \/pi / pi is the sound speed in the ionized gas. We note that due 
to the dynamic nature of this problem, the true solution should lie between 
the two curves (|68|) and ([69]) . since both radiation and gas pressure play a role 
throughout the dynamics, and neither T% or T e are in fact constant behind or 
beyond the I-front. 

We use the following problem parameters: the initial gas temperature is set 
to T = 10 2 K; the initial radiation energy density is E° = 1CP 20 erg cm~ 3 ; 
the hydrogen is initially fully neutral (i.e. HII/H=0); the spatial domain is a 15 
kpc box. We run for a simulation time of 1 Gyr, which is not long enough to 
reach final the final equilibrium rf, but well past the transition from R-type to 
D-typc. The implicit solver parameters are identical to those in £)4.51 but with 
a linear solver parameter St = 10 _9 ||/(C/fc)||. 

Results from these tests are shown in Figure [12] which plot the computed 
I-front position and neutral fractions for various spatial meshes, along with the 
"error" in these quantities. Since this test problem does not have true analytical 
solutions, we compute the "error" as the deviation in each solution from the 
most-refined 128 3 mesh solution. 

We also use this problem to examine the effect of our operator split solution 
strategy on the temporal accuracy of the solver. In Figure [T3] we plot the 
spherically-averaged temperature profile for a 128 3 spatial grid at 175 Myr, 
and the associated relative errors found through varying the time step size. The 
error plot has been zoomed in around the heated region and front. We note that 
although both PPM and the described implicit sub-solver are both up to second- 
order accurate, the splitting reduces the resulting accuracy to slightly better 
than first-order in time. We also note that the adaptive time-stepping strategy 
from Section 13.2.31 results in average time steps of 6.0e-4, 5.7e-4, 4.0e-4 and 
2.7e-4 for the 16 3 through 128 3 grids, respectively; whereas the hydrodynamic 
CFL-limited time steps for these same grids are 3.1e-3, 1.4e-3, 6.0e-4 and 2.8e- 
4. Hence the radiation and ionization dynamics drive the system for coarser 
meshes, while at finer mesh sizes the hydrodynamic CFL condition begins to 
dominate. 
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Figure 12: Convergence of the hydro dynamic ionization test with mesh refinement. Left: 
overlay of computed I-front position (top) and error (bottom) for varying mesh sizes. Right: 
overlay of computed neutral fraction (top) and error for varying grids. 
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Figure 13: Convergence of temperature profile as the time step size is refined. Left: the 
spherically-averaged temperature profile at 175 Myr. Right: relative errors in temperature for 
varying time step sizes. Average relative errors were 2.6e-4, 1.0e-4, 4.9e-5, 2.4e-5, 8.2e-6, and 
2.7e-6, for the coarsest to finest time steps, respectively. 
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Mesh 


Processors 


Time Steps 


Run Time 


Newton Its 


CG Its 


MG V-cycles 


64 3 


1 


266 


1694.38 


322 


914 


2991 


128 3 


8 


265 


2299.60 


274 


799 


2575 


256 3 


64 


265 


2456.58 


268 


787 


2524 


512 3 


512 


264 


2594.50 


265 


780 


2510 


1024 3 


4096 


264 


2707.30 


265 


780 


2510 



Tabic 2: Cosmological Weak Scaling Statistics. 



4-8. Weak Scaling 

As described in our introduction and throughout the description of our nu- 
merical methods, a key goal in introducting a fully implicit solution mechanism 
for the stiff components in radiation, hydrodynamics and chemical ionization 
simulations is the eventual scalability of such a solver to very large problem 
sizes. We therefore investigate the weak scaling of the implicit solver on the 
cosmological radiative ionization problem from section l4~6l For these tests, we 
emulate the setup from the q = 0.5 and Zj = 4 test, but here we place an 
ionizing source in the center of each processor's subgrid. Moreover, for these 
weak scaling tests we increase the domain size and mesh size proportionately to 
the number of processors, where each processor owns a 64 3 grid and an initial 
subgrid box size of 80 kpc. We then run problems that scale up from 1 to 4096 
processors, resulting in spatial grids ranging from 64 3 to 1024 3 . Moreover, since 
we are investigating the scaling properties of the numerical methods, we shorten 
the simulation time to evolve from Zj = 4 to z = 3 in order to conserve on su- 
percomputer resources, while retaining the portion of the simulation with the 
most rapidly-evolving dynamics. All runs were performed on the NSF Kraken 
machine (using 2 cores/node). 

We show the runtimes associated with these tests in Figure [TJ1 and provide 
detailed statistics from each run in Table [2j We note that on this architec- 
ture, the solver demonstrates near-perfect scalability, with modest increases in 
runtime for parallel versus serial runs, and only marginal increases in solution 
time as the parallelism is increased. The reason for this is the near-constant 
number of iterations required by the nonlinear Newton solver, the outer CG 
linear solver, and the inner multigrid prcconditioner. Therefore the increase in 
run time may be directly attributed to the ideal O(logp) increase in runtimes 
typical of multigrid methods, allowing near-optimal scalability to the limits of 
modern supercomputer resources. 

5. Conclusions 

We have described an implicit formulation for coupling cosmological radi- 
ation transport, chemical ionization and gas energy feedback within Enzo hy- 
drodynamics simulations. The formulation is based on an operator-splitting be- 
tween the non-stiff hydrodynamics and stiff radiation- ionization-energy feedback 
physical processes, in which the stiff processes are solved within a fully-implicit 
Newton-Schur-Krylov-Multigrid framework. 
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Weak Scaling, Cosmological Ionization Test (Kraken) 
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Figure 14: Weak scaling results for the cosmological HH-rcgion expansion test. 

Through numerous tests, we have demonstrated that this solver is accurate 
and stable, allowing simulations of a wide variety of physical environments from 
the laboratory scale to the astronomical and even cosmological scales. More- 
over, through the choice of numerical methods that form the implicit solver, 
it demonstrates ideal scalability for such coupled physics simulations. In addi- 
tion, this implicit formulation is highly extensible, and may easily be adjusted 
to allow new physical processes such as magnetic fields, multi-frequency radi- 
ation transfer, and additional chemical species. Finally, we are in the process 
of extending this approach to allow for adaptive spatial discretizations (AMR), 
which should require adjustments to only the inner multigrid linear solver. 
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